%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 09.22.23. Device 1 data (Half-Wave Dipolar Metal-Semiconductor Laser, Fig. 3b)
%%% Sangyeon (Fred) Cho, fredericksy.cho@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
lambda=open('lambda_device1.mat').lambda;
spectrum=open('spectrum_device1.mat').spectrum;
pump=open('pump_device1.mat').pump;
x0=300;
y0=300;
width=200*0.8;
height=300*0.8;
n = input('Enter a number: ');
switch n
    case 1 %Graph 2D
for i=1:21
figure(1)
plot(lambda,spectrum(i,:))
xlim([1050 1650])
xlabel('Wavelength (nm)')
ylabel('Ouptut (Photons)')
set(gcf,'position',[x0,y0,width,height])
set(gca,'TickDir','out'); 
hold on
end
    
case 2 %Smoothing 
for i=1:21
spectrum_s(i,:)=smooth(spectrum(i,:),'rlowess',100);
figure(1)
plot(lambda,spectrum_s)
xlim([1050 1650])
xlabel('Wavelength (nm)')
ylabel('Ouptut (Photons)')
set(gcf,'position',[x0,y0,width,height])
set(gca,'TickDir','out'); 
hold on
end

[X, Y] = meshgrid( pump,lambda);
figure (2)
plot3(X, Y, spectrum_s);
ylim([1100 1600])
set (gca,'Ydir','reverse')
set(gca, 'colormap', jet)
set(gca,'TickDir','out'); 
ylabel('Wavelength (nm)')
zlabel('Ouptut (Photons)')
xlabel('Pump(mJ/cm2)')

case 3 %Input-Output
for i=1:21
m=spectrum(i,:)-spectrum(5,:).*mean(spectrum(i,1300:1500))./mean(spectrum(5,1300:1500));
mode(i,:)=m;
figure(1)
plot(lambda, mode(i,:))
xlim([1050 1650])
xlabel('Wavelength (nm)')
ylabel('Ouptut (Photons)')
set(gcf,'position',[x0,y0,width,height])
set(gca,'TickDir','out'); 
hold on
out(i)=sum(mode(i,650:950),2);
end
figure(2)
plot(pump.*10^3, out,'o')
xlabel('Pump (mJ/cm2)')
ylabel('Ouptut (Photons)')

    case 4  %Linewidth
        for i=1:21
m=spectrum(i,:)-spectrum(5,:).*mean(spectrum(i,1300:1500))./mean(spectrum(5,1300:1500));
mode(i,:)=m;
        end
for i=1:7 
F=fit(lambda(470:1680).',squeeze(spectrum(i,470:1680)).','gauss1');
fwhm(i)=F.c1.*sqrt(2.*log(2));
end
for i=8:21
F=fit(lambda(650:850).',squeeze(mode(i,650:850)).','gauss1');
fwhm(i)=F.c1.*sqrt(2.*log(2));
end
figure(1)
plot(pump(8:21).*10^3, fwhm(8:21),'o')
xlabel('Pump (mJ/cm2)')
ylabel('Linewidth (nm)')
end

